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Abstract. The continued-fraction method to solve classical Fokker-Planck 
equations has been adapted to tackle quantum master equations of the Caldeira— 
Leggett type. This can be done taking advantage of the phase-space (Wigner) 
representation of the quantum density matrix. The approach differs from those 
in which some continued-fraction expression is found for a certain quantity, in 
that the full solution of the master equation is obtained by continued-fraction 
methods. This allows to study in detail the effects of the environment (fluctuations 
and dissipation) on several classes of nonlinear quantum systems. We apply 
the method to the canonical problem of quantum Brownian motion in periodic 
potentials both for cosine and ratchet potentials (lacking inversion symmetry). 



PACS numbers: 05.40.-a, 03.65.Yz, 05.60.-k 

1. Introduction 

The phase-space formulation of quantum mechanics has received a renewed attention 
in the last decades 12 0] because it allows to employ notions and tools of classical 
physics in the quantum realm. The central object in this approach is Wigner 's function 

Wi^'P) = ^ J dye-'Py/^gix+^y,x-^y) (1) 

which is merely a phase-space {x,p) representation of the density matrix g{x,x') — 
{x\g\x'). In a closed system, the dynamical equation for the Wigner function is 
equivalent to Schrddinger's (or von Neumann) equation while, in the ?i — > limit, 
it recovers Liouville 's equation for the phase-space distribution in classical mechanics. 

Another remarkable property of the Wignerian representation is the average 
property. The expectation value of a quantum operator A is obtained from the 
corresponding classical variable A{x,p) (related with A by Weyl's rule ^ Ch. 1]) 
by a prescription analogous to the classical one, namely 

(A) =Tr(^i) = J dxdpW{x,p)A{x,p) . (2) 

Thus, in spite of not being necessarily positive, the Wigner function can be considered 
as a sort of quantum-mechanical "distribution" . Besides, the marginal distributions 
of W give the true quantum probabilities of x, as P{x) = {x\g\x) — J dpW{x,p), or 
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of p, as P{p) = {p\q\p) — ^ dxW{x,p). Therefore, the Wigner formahsm provides a 
natural quantum-classical connection. This has been specially exploited in the search 
for quantum analogues of various classical effects (e.g., chaos) or in the problem of the 
fuzzy borders between the quantum and classical worlds [S]. 

In contact with the surrounding medium, the (open) system experiences 
dissipation, fluctuations, and decoherence O 13 El- Typically one is interested in 
a system consisting of a few relevant degrees of freedom coupled to its environment, 
or bath, which has a very large number of degrees of freedom (e.g., phonons, photons, 
nuclear spins, etc.). The system is not necessarily microscopic, but it can be a 
mesoscopic system described by some collective variables which under appropriate 
conditions can behave quantum mechanically 0^]. Because of the genericity of these 
situations and the fundamental issues raised (e.g., approach to thermal equilibrium, 
dissipation in tunnelling and coherence, measurement problem), the study of quantum 
dissipative systems is of interest in several areas of physics and chemistry 

The dynamics of an open system can in many cases be formulated in terms 
of a quantum master equation for the (reduced) density matrix. In the Wigner 
representation this can be compactly written as dtW = CW, with £ a certain 
evolution operator. Taking the classical limit, the master equation for a particle 
of mass M in a potential V{x) goes over the Klein-Kramers equation |E1 ^| 

dtW=[- (p/M) d^ + V'dp + j dp {p + MkBTdp)] W (3) 

which is simply a Fokker-Planck equation in phase space. The first two terms (Poisson 
bracket) generate the classical reversible evolution (Liouville equation). The last 
term ( "collision operator" ) accounts for irreversible effects due to the coupling to 
the environment (dissipation and fluctuations). The damping parameter 7 measures 
the strengh of the coupling to the bath, which is at temperature T. It is worth 
recalling that phase-space dynamics also includes problems reducible to "mechanical" 
analogues: Josephson junctions, certain electrical circuits, chemical reactions, etc. 

A suitable non-perturbative technique to solve classical Fokker-Planck equations 
of systems with few degrees of freedom is the continued- fraction method |14| . This is 
a special case of the expansion into complete sets (Grad's) method to solve kinetic 
equations in statistical mechanics I15j . In brief, the distribution is expanded into a 
series of orthogonal polynomials and, as a result, the kinetic equation is transformed 
into an infinite set of coupled equations for the expansion coefficients Ci . Approximate 
solutions are obtained by truncating the hierarchy of equations at various levels. To 
get manageable expressions, however, the truncation needs to be performed at a low 
level. In the continued-fraction method, instead of truncating directly, one seeks for a 
basis in which the range of index coupling is short (ideally, the equation for Ci involves 
Cj_i, Ci, and C^+i). Then, the (differential) recurrence relations among the Ci can be 
solved by iterating a simple algorithm, the structure of which is like that of a continued 
fraction | |Appendix A| ). For the Klein-Kramers equation (PJ the tridiagonal chain of 
coupled equations for the Ci is called the Brinkman hierarchy .16i) . It has been solved 
by continued-fraction methods to study classical particles subjected to dissipation and 
fluctuations T4'. The method has been also extended to rotational Brownian motion 
problems involving classical spins and dipoles (vd. |17j and references therein). 

To deal with quantum dissipative systems, however, is a more delicate task [7j. To 
begin with, phenomenological quantisation of dissipative systems poses fundamental 
problems (e.g., with the uncertainty and superposition principles). A rigorous route 
is to model the bath in a simple way, quantise it together with the system, and 
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eventually trace over the environmental variables. However, the resulting theoretical 
frameworks are technically involved. Quantum master equations, except for simple 
cases, are difficult to solve. An alternative is provided by quantum state diffusion 
methods (vd. 2HI)i where stochastic evolution equations for state vectors in Hilbert 
space are introduced as computational tools. Nevertheless, its implementation seems 
to be restricted to systems with discrete spectrum (e.g., oscillators, 2-state systems, 
etc.). Similarly, quantum Langevin equations are of limited use beyond nearly 
harmonic systems. For an arbitrary system, there exist exact path-integral expressions 
for the evolution of the density matrix (involving the so-called influence functional that 
incorporates environmental effects) 17] . However, those expressions are difficult to 
evaluate, even numerically, because the propagating function is highly oscillatory, 
rendering numerical methods unstable at long times. Finally, quantum Monte Carlo 
simulations can in principle be used. Nevertheless, in spite of ongoing progress (vd. 
|18|'). they are computationally complex and suffer from the (dynamical) sign problem. 

This situation strongly motivates the development of alternative methods for 
quantum dissipative systems. Inspired on their suitability for classical systems, 
continued-fraction techniques have been developed for a number of problems. Shibata 
and co-workers j2()l |2] applied them to solve a c-nuniber quantum Fokker-Planck 
equation for a spin in a dissipative environment. Vogel and Risken employed 
continued-fraction methods to solve master equations in quantum nonlinear optics (vd. 
also [331). Recently "57 this method has been extended to study genuine phase-space 
problems within the Wigner formalism. A generalisation of the Brinkman hierarchy 
for quantum master equations of the Caldeira-Leggett type was presented. It was 
shown that the continued-fraction method for the classical problem can, in principle, 
be adapted to solve this hierarchy, yielding a promising technique to study several 
classes of nonlinear quantum systems subjected to environmental effects. 

In this article we first present the details necessary for the derivation of the 
quantum hierarchies and discuss their solution by continued fractions. The approach is 
then implemented for the problem of quantum Brownian motion in periodic potentials 
(a demanding problem with (partly) continuous spectrum). Both the cosine and simple 
ratchet potentials (lacking inversion symmetry) are considered. For the former a 
number of previous results are recovered (so testing the method) and extended, while 
the physical interpretations are revisited within the Wigner formalism. For particles 
in ratchet potentials we study the effects of finite damping (inertia) on the rectified 
velocities. The phenomenology is interpreted in terms of the interplay of thermal 
hopping, overbarrier wave reflection, and tunnelling. Results for non-equilibrium 
dynamics under oscillating forcing are also discussed, which show the combined effect 
of quantum phenomena, thermal activation, and dissipation in a nonlinear system. A 
number of technical issues are consigned to the appendices. 



2. Quantum master equations in phase space 



We start from the following generic form for the quantum master equation in the 
Wigner representation PSI 051 ITzl 05] : 



s— 1 ^ ^' 



W . (4) 



The first three terms, identifying Dpp — jMkBT, give the classical Klein-Kramers 
equation The mixed diffusion term DxpdIpW is heuristically related with the 
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colour of the quantum noise j^. The series of derivatives Sp^^^^'VF (Wigner-Moyal 
term) gives the quantum contribution to the unitary evolution of the closed system. 
Conditions under which this series can be truncated are sometimes discussed, to avoid 
dealing with an infinite-order partial differential equation. However, we shall fully 
keep this term because it can be specially important in nonlinear systems 

Conditions of validity for the master-equation description are discussed in |81ll!i2l 
1^ . Equations of the type of Q can be derived modelling the particle surroundings 
as a hath of oscillators representing the normal modes of the environment. However, 
a number of assumptions are typically involved like semiclassical or high-temperature 
bath (rendering the kernels of the influence functional local in time), weak system-bath 
coupling (Born-Markov approximations), etc. Then, one would heuristically expect 
equation being valid for small enough h'j/kBT. Notwithstanding this, the structure 
of the equation seems to be quite generic. Thus, a quantum master equation recently 
derived for strong coupling and valid for all h'-f/ksT, involves terms similar to 
those in (0J [with x-dependent coefflcients and renormalised V^(a;)]. Besides, the phase- 
space representation of the celebrated Lindblad master equation is obtained by simply 
adding to Q terms of the form dx{xW) and d^W Most of these terms can be 
readily incorporated in the treatment below (indeed, the term Dxpd^pW, absent in 
the original Caldeira-Leggett equation, is included to illustrate this). The same, in 
principle, would apply to possible extensions of the quantum master equation (0J. 



3. Preliminary manipulations of the master equation 

First it is convenient to introduce appropriate scaled units. These are based on 
a characteristic length xq (e.g., the distance between the minima in a double- well 
oscillator, the period in a periodic potential) and a characteristic energy Eq (e.g., a 
barrier height). Then one scales frequencies by uiq = {Eo/Mx'^y^'^ , time by ijJq^, forces 
by Fq = Eq/xo, the action by So = Eq/loo, etc. On the other hand, Dpp {= jMIcbT 
in a high-T bath) is handled as an effective temperature kBT^s = Dpp/jM and we 
introduce some convenient thermal rescalings via ut = (fcBTeff/Ma;g)i/2. Thus, in the 
equations p is scaled by Mxqujt, the potential and D^p appear divided by kBT^s, 7 
enters in the combination 77^ = 7/^T: time is multiplied by lot, while the thermal de 
Broglie wave length A^b = V(4MfcBTcg)'^/^ enters divided by Xq. 

Omitting any marks for the scaled quantities (but keeping in mind specially the 
thermal rescaling of t, V{x) and p), the master equation is simply written as 

dtW^ [-pdx+v dp+jTdp{p+dp)+Dxpdlp+j:Zi^^"^ y^^"^'Hx)4'''^'^]w (5) 

with the coefficient in the quantum sum given by 

^^^^ = i~irxl%/{2s + iy.. (6) 

Planck's constant h is introduced in terms of the characteristic action Sq = Eq/luq via 
the quantum parameter K (denoted K in j24|') 

h/So = 2TT/K, AdB=7r7T/a, a={'j/uJo)K. (7) 

The second relation gives AdB in terms of the Kondo parameter a, related by the third 
one with our K. The classical limit is approached as H/Sq 0, i.e., letting K 00. 

Let us now introduce a splitting of the evolution operator £ that will facilitate 
the calculation of the matrix elements required when expanding W{x,p) into complete 
sets. On inspecting equations ((SJ and ©, we see that extending the quantum sum 
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down to s = 0, we obtain the part V dp of the classical term. Thus, we can decompose 
C into the irreversible, kinetic, and potential parts: 

{-Cir = 7t dp (p + dp) 
/Ikin = -{p- D^pdp)do, . (8) 

This natural splitting has the advantage of dealing with V(x) (the most problem 
dependent part) in a unified way. Besides, we have grouped the term dpdx with the 
kinetic term pdx because they are structurally similar. 

4. Derivation of the quantum hierarchies 

In the expansion into complete sets approach to solve kinetic equations, the 
distribution is expanded in an orthonormal basis {ijjn} and the coupled equations 
for the coefficients C„ derived |15[ p. 175]. For the Klein-Kramers equation Q one 
uses Hermite functions of p as basis [211 Sec 4.4] 

Mp) = e-p'/^H„{p/V2) , r„ = [(27r)i/22"n!] . (9) 

One starts with the expansion in the momentum because, while the p dependence 
of the Hamiltonian is fully specified (kinetic energy p^/2M), the potential part V{x) 
depends on the problem. Thus, explicit manipulations can be done on the parts 
involving p, which are valid for any system. On the other hand, the Hermite basis has 
a number of advantages ; not the least is the handling of the derivatives dp in the 
dynamical equation by means of the associated creation and annihilation operators 
b = dp + ^p and = —dp + ^p (this will prove very convenient in the quantum case). 

For the Klein-Kramers equation the resulting equations for the expansion 
coefficients C„ are called the Brinkman hierarchy 1 141 116| . This plays an important 
role in classical systems, both for analytical treatments — derivations of 1 /j expansions, 
etc. — and to develope efficient non-perturbative numerical methods. In fact, it has 
the structure of a 3-term recurrence relation which, after expansion in the position 
basis, can be solved by continued fractions. In what follows we shall proceed in the 
quantum case following as closely as possible the steps of the classical limit. 

4.I. Expansion in the momentum basis 
We start expanding the Wigner function as 

W{x,p) = wJ2Cn{x)i^n{p), "' = 7^^^"*^'^ (1"^) 

where we have extracted the Boltzmann-type factor w. This involves the auxiliary 
parameter < rj < 1/2 and potential ^{x), which is frequently chosen proportional to 
V{x), i.e., $ = eV. The results should not depend on the 77 or e used, but these may 
(i) put parts of the evolution operator C in self-adjoint form and/or (ii) improve the 
stability and convergence of the eventual numerical implementation. 

From the orthonormality of the ipnip), the expansion "coefficients" can be written 
as C„ — Jdpipn dm'i'm — J <ip i'n (W/w). Then, differentiating with respect to t 
and using the evolution equation dtW — C W, we get the dynamical equations for C„ 

dtCn = y^^QnmCni , Qnrn ^ dpipn^'fpm, C=W^^Cw. (11) 

m •' 
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m ; 



To get the matrix elements Qmm which are stiU operators on the x-dependence of 
one needs the w~^{ ■ )w transform of C . This can be spht as £ = C\r + /^kin + 'Cv, 
with £ir, £kin, and >Cv given by ((HJ- The transformation of these operators is done in 
[Appendix B| taking advantage of resuhs for normal ordering of h and 6"*" (which we 
use to express the dp and p's in £). Then, with the C„, £kin, and £v obtained we 
compute their matrix elements between the ?/'ti, getting the Qnm- The details of this 
part of the calculation are given in [Appendix C| 

Introducing the results obtained for Qnm into dtCn ~ X)™ QnmCm one gets the 
following quantum (Brinkman) hierarchy 

[(«-l)/2] 



=0 



+ ^/^.D-Cn-l+lnCn + ^frlTl D+Cn+1 (12) 

oo 

+ v/(" + l)(n + 2) 7+ C„+2 + [r:^+ ^^''+'^] C^n+(2.+i) • 

The auxihary (damping) parameters introduced read 

7± = 7T?7±(1 - ?7±) , 7„ = 7T [2n(?7 - <t) - (13) 
which involve the following ?7-related parameters (note the sign exchange) 

V±^VT^, <^^V-V+- (14) 
In equation (|12|l . the operators on the x dependences are 

= d± (9, - , d± = 1 + r;±i?,p , A = aXl^dl , (15) 
T'n^ = vl'^'iil:^ c-^/^iFi ( - m, 2s + 2 ; A) , = k^^) v^^^iA^ ■ (16) 

Here m — n~ (2s + 1) , i-Fi (a, c ; z) is the confluent hypergeometric (Kummer) function 
[ 371 Ch. 13.6], and k'*' is given by equation ®. In is important to note that the 
act only on V{x) not on the Cm(x) ( [Appendix B[ |. 

The quantum hierarchy p2|l is equivalent to the Caldeira-Leggett master equation 
(@J. Previously, several hierarchies had been derived ' 381 1391 1^ 141 j . but the treatment 
was for closed (Hamiltonian) systems J8, 39, 40 or involved unsuitable bases like {p"} 
[ 381 14l)[ |41 [ . In contrast, our hierarchy of Hermitian moments constitutes a direct 
quantum generalisation of the Brinkman hierarchy. To verify this, note that in the 
classical case only the s = terms survive in the sums of lO ( [Appendix Q . Then, 
setting the usual rj = 1/2, equation H12|l reduces to 



dtCn 



{D- + V') Cn-i + n7T Cn + D+Cn+1 (17) 



which is indeed the celebrated hierarchy associated to the Klein-Kramers equation. 

For potentials obeying V^^^ = 0, Vs > S, the quantum hierarchy (|12|) is a 
recurrence relation with coupling to a finite number of terms. Examples include 
the harmonic oscillator V = ^Mlu^x^ and the double- well (Duffing) oscillator V = 
— ^ax^ + jbx'^. Finite-coupling recurrences are always reducible to 3-term vector 
recurrences ( [Appendix A[ ) and solvable by continued fractions, as in the classical case. 
However, for non-polynomial V{x) this approach is prevented by the infinite coupling 
range in the index n, consequence of the infinite series of p-derivatives T/^^^+i) dj^'^^^^ 
in the Wigner-Moyal term. In what follows we shall exploit the expansion in the 
position basis to show how this problem can be circumvented. 
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4..2. Expansion in the position basis 

Recall that the coefficients C„ are still functions of x. Inserting their expansion in an 
orthonormal basis {ua{x)} into the quantum hierarchy we get the generic form 

with matrix elements {B)ap — Jdx u'^ B(x, dx) up. To put H12I) into the above form 
we just need to (i) substitute the operators (and x-dependent coefficients) by their 
matrix elements {B)ap, (ii) attach position superscripts to the C„, and, finally, (iii) 
sum over them. On so doing we get the expanded quantum hierarchy 

- -EE [rr ^^"+^'] c^l(..+r) + yF3T>^ E 7- s^, ct. 

s>0 P 13 

C^ + V^E^o/jCf+i (19) 

13 13 13 

+ Vin + l)(n + 2) ^ 7+ S^p + E E Ki2s+i)y^''^'^].p • 

Considering that = d±{dx — ^') and $ = eV, we see that to get all matrix elements 
involved we just need those of dx and of y(2s+i)^ namely 

{d.)^p = jdxuldxup, V^f^'^ = JdxulV^'^+'Kf, (20) 

because the only increase the order of the potential derivatives. 

The expanded hierarchy has the form of a system of ordinary differential 
equations, — J2m0Qnm^m- If ^ A are some large truncation indices for 
the bases {tpnip)} and {ua{x)}, we have a problem with (iV x A) equations. One may 
be tempted to try direct Runge-Kutta integration or numerical diagonalization of the 
associated {N x A) x {N x ^4) matrix. However, the dimensions involved are typically 
very large, so that it is worth to pursue a continued-fraction treatment. | 

As mentioned before the p-recurrence has a finite coupling range for polynomial 
potentials. For other potentials, with an appropriate choice of the basis {ua{x)}, the 
matrix elements (|20f) may vanish when the second index /3 lays at a certain "distance" 
of the first a. Then, the expanded hierarchy would be a recurrence relation in the 
position indices with a finite coupling range, also amenable for a continued-fraction 
treatment. To proceed in these cases, we just need to transform the two-index extended 
hierarchies into ordinary recurrence problems (i.e., with one index). 

5. Matrix quantum hierarchies with one-index recurrences 

In order to apply the continued-fraction method when finite coupling range in one 
of the indices is attained, we need to convert first the two-index recurrences into 
ordinary recursions. In the resulting one-index recurrences the coefficients will be 
matrices acting on appropriate vectors formed with the C". 

I Dynamical equations for the Wigner function have anyway been tackled by purely numerical 
methods, like (pseudo-) spectral methods for partial differential equations (see 1421 1431 and references 
therein), or grid discretization usually followed by Runge-Kutta integration 1441 . 
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The generic expanded form H18I) can be converted into an one-index recurrence relation 
in the momentum index introducing the following vectors and matrices 

/ (Qjim)_^ _^ • • • (Qrmi)_^^ ^ 

: ■.. : (21) 



dt 



I A, -A 



('3"'") A, A / 



where iQnm)ap = /da;u; 
with all the C" with a 



Up. That is, for a fixed n one forms the column vector 
= —A, . . . , A (this index may take positive and negative values). 
Similarly, for given (n, m) we build the matrix Q^m = {Qnm)ap with all indices a, (3. 

Before giving explicit expressions for the matrix coefficients we introduce the 
following notation: = Q„,„±„ Q±± = Q„,„±2, Q± = Q„,„±i, and Q„ = Q„,„. 

Then, the matrix quantum hierarchy H21|l can be written as 

d 

d^^ 

"„c„ + Q+c„+i (22) 



'Cn+2 + J2 



s>l^n 



ll+(2s+l). 



'n+(2s+l) 



We have accounted for Qt^^'^^ = 0, Vs > 2 [vd. equation (fT^ ] and incorporated the 
s = terms into the central c„±i ones. The matrix coefficients explicitly read 

/^-(2s+l)\ _ 



'q/3 



■ )q 



'q/3 



/3 

-"''') a/3 



/3 



- l)n 7_ 6a(3 

L « J a/3 
— 7" <5a/3 

[r:;:i^']a/3 



(23) 



a/3 



- V(" + !)(" + 2) 7+ '5a/3 



'a/3 



- r 



V^(2s+l)l 



ri+(2s+l) ' J a/3 

For finite coupling range in n (e.g., in polynomial potentials), the hierarchy (|22l) is 
a recurrence relation to which continued-fraction methods can now be applied. This 
amounts to convert the underlying large {N x A) x (N x A) problem into N problems 
with associated Ax A matrices [or A {2A+ 1)]. As the recursion "coefficients" are 
matrices, one would need matrix continued-fraction methods ( [Appendix A| ). 



5.2. Matrix quantum hierarchy: x-recurrence 

When the p-recurrence has an infinite coupling range (in periodic potentials, the 
Morse potential, etc.), we can convert the 2-index hierarchy (|18|) into matrix form, 
but introducing recurrences in the position index a 



At 







la!3 



N 



( {Q 



Q/3 



a/3 



(Q 



NN I 



(24) 



a/3 



That is, (Qa/3)nm = {Qnm)af}, with Baft = jdxu*^Bup. Thus, for a fixed a one 
constructs the column vector with the for all n = 0, . . . , TV, and for given (a, /3), 
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the matrix Q^m 
equation (|23|) for 



(Onm)a/3 with all n, m. 



we can construct 



Ja/3)„,„_1 = 
Ja/3)„„ = 

^"^)n,n+(2s+l) ^ 
all other diagonals being zero (C^a/3)n,n=p2 



Comparing the matrix in ()24|l with 
ap diagonal by diagonal 

- ^(n - l)n 7_ (5a/3 



- v^(n+ l)(n + 2) 7+ (5„/3 



(25) 



0, Vs > 2. 

To recognise better the structure of these matrices, we decompose them into the 



"free" and potential contributions 



lap 



tp 



tap- 



For 



including the kinetic 



and irreversible terms, we have the pentadiagonal structure (tridiagonal for = 1/2) 
/ 



tap 



iQ^aP 




s/l^'^+SaP 











^D-p 


ll^aP 












^T^-ySap 


^D-p 


j2SaP 


^Kp ^ 


y/S^-y+Sap 








'/2^^-5ap 


^D-p 


IZ^aP 




y4^7+5a/3 








s/i^-ySap 




JiSaP 

















Ib^ap 



The part due to V{x), on the other hand, has the following alternate dense structure, 
reflecting the odd powers of dp in the Wigner-Moyal term: 



lap 












aP 



[T\+V(^^^p [vl^v^'Aap 











[r: 



-v%p [rr^Vp 

[Tl^-V^^)U [Tl-V%p 

[r^'"y(3)]„/5 

Wi-v^-'n.^a 



[Tl'^V^'\p 
y'\aa 

[r°'+V']a/3 



With an appropriate choice of the basis {^^(a;)} the matrix elements {B)ap may couple 
only a finite number of basis functions. Then, the recurrence (|24l) would be solvable 
by continued fractions when the rt-recurrence is infinite and this approach is seemingly 
precluded. Even for finite coupling in n, the a-recurrence may be preferable in certain 
situations (e.g., at low T), as indicated by Risken for a cosine potential in the classical 
limit ^1 Sec. 11.5.6]). In these cases, one has been able to reduce the underlying 
large (A'' x A) x {N x A) problem to A problems with associated N x N matrices. 
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6. Density matrix, observables, and marginal distributions 

Once the c„ = QnmCm or Cq, = J^p Qafscp are solved, we can reconstruct the 
Wigner function from its expansion coefficients C" in the double basis {ua{x) ipnip)} 

W{x,p)=w{x,p)Y,C:uo.{x)Mp) (26) 

obtaining the full solution of the quantum master equation. If preferred, one can switch 
to the familiar density matrix g by the inverse of transformation 

g{x,x') = / dpe'P'-'^~'''^^^W{^{x + x'),p) . (27) 



We can get any observable inserting the above W{x,p) in equation for the averages. 
Nevertheless, common observables can many times be extracted directly from the C". 

Let us illustrate this with transport observables, which are characterised by the 
averages {p^}- For an arbitrary function of p only, we get from (|26|l 

(/) - / dxdpWix,p)fip) ^Y.^nf dxe-^^-^u^ix) [ dproe-'^P"/'fip)Mp) (28) 



where we have used w = r^e ''P^/^e *, with rg — l/(27r)^/* [equations ^ and (|10|l ]. 
Next we introduce the auxiliary integrals (for their calculation vd. [Appendix D| ) 

dxe-*(-)7.„(a;) , ifW = /dpro e-^^'/^ / Vn(p) • (29) 



Then, the moments (case / = p^) can be written in terms of the expansion coefficients 
as (p^) = J2n ( ^a) ■ Taking into account that _ff„ has the parity of n and 
carrying this to Kn^ , we get for the first two 



The first moment is the current and the second characterises fluctuations around it. 
The zeroth moment corresponds to the normalisation of W, and imposses on the 

the condition (p") - E„i^2"HEa t;2"„^a) - 1- 

Finally, the quantum probabilities of x and p, given by the marginal distributions 

P{x) — JdpW{x,p) and P{p) — /da; T4^(x,p), can be expressed in a similar fashion 
P(x) = e-^(-) E ^"(^) ( E ^" ^n^) (31) 

a n 

P{p) = roe-'^^/^E^«(^^) (E^" ^0 • (32) 



7. Periodic potentials: matrix elements and limit cases 

Henceforth we shall apply the method described to incorporate fluctuations and 
dissipation in the problem of quantum transport in periodic potentials. The simplest 
model consists of a particle evolving in a cosine potential subjected to external forces 
(tilted "washboard" potential), which is also a paradigm of quantum Brownian motion 
[ 451 146L IT7| . Others include potentials lacking inversion symmetry (ratchets), which 
have been used to model rectification of current and directional motion in several 
systems 02] . These problems are demanding because periodic potentials have (partly) 
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continuous spectra |49| . rendering inappropriate methods devised for systems with 
discrete levels. In this section we compute the matrix coefficients of our recurrences 
for arbitrary periodic V{x). With them we can implement the continued- fraction 
method to solve the quantum master equation, which we first check by regaining the 
classical dissipative and quantum Hamiltonian limits. 



7.1. Matrix elements 

In periodic potentials the p-recurrence cannot be used because of its infinite coupling 
range (the derivatives of V{x) neither vanish, nor decrease, as they are essentially the 
potential itself). However, one can guess that the x-recurrence may have a coupling 
range related with the number of harmonics in V{x). To exploit this we introduce 
plane waves as basis functions and the Fourier expansion of the potential derivative: 

u^{x)^^, y'(x)=5]F>--. (33) 



To preserve the periodicity of W , one extracts the external force F from V{x) including 
it in a generalised — d±[dx ~ $') — rj±F] cf. equation (|15|l . Then, the auxiliary 
potential <& is set proportional to the periodic part $ = eV . 

To compute the matrix elements in Qq,^ , we only need to do the integrals H2Q(I with 
plane waves for Ua. Using dxUa = ia e'"^/-y27r, the integral J^^^'^da; e'*^'^~"'^ = 27^5^(3, 
and the expansion V^'^'' = vy\^^^ , we find 

{dxU = ia 5^0 , [V<^'^ {x)] „^ = . (34) 

From these elements and Vq^'^^^'' = {—lYg'^^V' we get 



g 



(35) 



D^p = (}d±a--q±F) S^p - d±e V'^_p 

with q — a — (3,m — n~ (2s + 1) and [cf. equation Hlt)|) ] 

G'niz) = |4^)(9)|e-^/Si^i(-m,2,s + 2;z) , |4^)(g)| = ^^^M . (36) 

(2s + Ij! V m\ 

That is, Kn\q) is the coefficient in (|16f) with AdB — * 'Z'^dB- Inserting equation (|35|l 
into the general matrices Qap of section Owe explicitly get the matrix coefficients of 
our recurrences ( [Appendix E| ). For the cosine potential 

l/(x) = -Vbcosx (37) 

the coupling range is 1, because = for |q| > 1 and we have a 3-term recurrence 
Ca = Qa,Q-iCQ--i + Qa.aCa + Qa.a+iCa+1- For a 2-harmonic ratchet potential 

V{x)^-Vo[ sin x + (r/2) sin(2x)] (38) 

the range is 2 (5-term recurrence) because = for |g| = \a — f3\ > 2, and so on. 

As we restricted < r] < 1/2, we have a = ??-77+ — — 1/4 < 0, so that 
z = — crg^A^B > 0. Thus, the argument of is positive and exp(— z/2) can act as a 
regularisation factor [an advantage of allowing 77 ^ 1/2 in (|10l) ] . This factor reduces the 
weight of the off-diagonal terms inside Qa/3, enhancing the numerical stability when 
going into the deep quantum regime. Then we typically use r] ~ 0-0.05, while for 
classical calculations 77 = 1/2 performs better. Concerning the auxiliary potential $, 
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Figure 1. Wigner functions in the classical limit [K = 1000) for a particle in 
a tilted sinusoidal potential (two periods are displayed along x). Here 7 = 0.05, 
and T = 1, while F = 0, 0.075, 0.15, and 0.2 (left to right, top to bottom). 

for the problems considered below we use £ = 0. Convergence is achieved with N ^ 100 
Hermite functions and ^ ~ 25 — 50 plane waves [then TV X (2A + 1) ~ 10^]. Finally, 
a word on the scaled quantities (section 121) . To get period 27r, we have set xo — a/2iT 
with a the original period. The characteristic energy is given by the potential 
amplitude Eq = Vq] then k^T is scaled by Vq and the other characteristic quantities 
are: action, Sq — Xq\/AIVo, frequency, Uq = Vq/Mx^, and force, Fq — Vq/xq. 

7.2. Classical limit 

Before going into the quantum regime, let us check that for small enough h/ So = 2tt/K 
we recover the classical results. Recall that in the deterministic limit a particle in a 
periodic potential has two critical forces Fi (retrapping force) and F3 (force at which 
the barrier disappears). For F < Fi the attractors of the dynamical system are 
solutions "locked" around the potential minima, whereas for F > Eg, the "running" 
solution is an attractor globally stable. In the range Fi < F < F3 the system exhibits 
bistability. Finally, Fi decreases with 7 [for the cosine potential Fi ^ (4/7r)VJ3 7], 
whereas Fi tends to F3 as the damping increases, narrowing the bistability range. 

Thermal fluctuations are incorporated by the Klein-Kramers equation, which 
for classical particles in periodic potentials was solved by continued fractions in 
^^innilSE2| (using the p- recurrence). Here we solve the Caldeira-Leggett quantum 
master equation using a large K. As shown before, when V{x) = —Vq cos(x) we 
have a 3-term recurrence which, for the static response, can be written as 

QaC„_i+Q„c„ + Q+c„+i = (39) 
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Figure 2. Energy bands for various K [reduced Kondo parameter J7J]. Tlie 
potential profile is plotted to show the number of bands below the barrier. Left 
panels: cosine potential for K = I (as in figures 1^ and El . with the free-particle 
parabolic relation E oc q-^, and for K = 200 (figure |7J, together with the first 
levels of the associated anharmonic oscillator [dots; equation I50i ] . Right panels: 
ratchet potential jSU with r = 0.44 for K = 10 and K = 20, as in figure l6l 



0.5 



with 



and 



^Q,a-fi- Solving it with the continued- 



fraction method we reobtain Risken's impressive graphs for the classical distribution 
W{x,p) |14[ Ch. 11.5]. Some of them are displayed in figure^lfor weak damping and 
various external forces (here Fi ~ 0.064 and F3 = 1). At low F the distribution, 
always periodic along a:, has maxima at the potential minima, while the profile in 
momenta is a Maxwellian envelope ~ exp(— As F is increased we see the 
evolution from these locked solutions ( ( p ) ~ 0) to bistability between the locked and 
running solutions, and eventually purely running states W ~ exp[— (p — F/j)'^/2] are 
obtained at large forces. 



7.3. Hamiltonian quantum case 

7.3.1. Bands and quantum parameters. Having checked the retrieval of the classical 
results, we now enter the quantum regime by decreasing K ~ Sq/H. Firstly, it 
will be useful to get a more quantitative feeling of the relation of our parameter 
K and "how quantal is the system" . To this end we study the bands of the 
corresponding Hamiltonian problem (without force). Inserting the Bloch function 
^mq{x) oc exp(iqx)C/mq(a;) in the Schrodinger equation, we arrive at the following 
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Figure 3. Wigner functions in tiie absence of dissipation for a particle in a 
sinusoidal potential (two periods are displayed). The reduced Kondo parameter 
is A" = 1 (left) and 2 (right). The "islands" correspond to zones of negative W . 



{iq + d^f + V{x) U{x)=EU{x) , -i < g< i . (40) 



eigenvalue problem (in scaled units; section |2Jl 
- l/27rN" 

Expanding U{x) in plane waves (|33|l and truncating \a\ at some large A one gets a 
{2A + 1) X {2A + 1) matrix which is diagonalised numerically. Various bands E{q) 
pertaining to problems discussed below are displayed in figure |21 We see that the rule 
of thumb is that K/2 is of the order of the number of bands laying below the barrier. 
Finally, the quantum parameter used by Kandemir is (K/tt)'^, while Chen and 
Lebowitz employed flq = n/ K |54[ I55| . Besides, for a Josephson junction problem we 
have K = 7r(i?j/2i?c)^^^, with Ej and Ec the Josephson and Coulomb energies. 



7.3.2. Wignerian representation. Let us continue for a while in the Hamiltonian 
case to familiarise ourselves with the structure of the Schrodinger eigenstates in a 
phase-space representation. Kandemir obtained analytical approximations for 
the eigenfunctions of a particle in a cosine potential in terms of Mathieu functions. 
They were represented through the Wigner function associated to the density matrix 
of pure states g{x,x') — '^{x)'^* {x'). (In [SS] various Dirac deltas occurring were 
regularised with Gaussians S{p) ~ (k/Tr) exp(— fc^p^) and in the plots k = 1 was used; 
we attain the same setting T = 1 since our p is thermally rescaled.) 

Using a very small damping (7 — 10^^) in the Caldeira-Leggett master equation 
we reobtain Kandemir's Hamiltonian results for the ground state (figure 13)1. We see 
that at low K the state is extended or delocalised, becoming more localised at the 
potential wells when K (x Sq/K is increased. According to the regions of negative 
W are associated to Brillouin zone boundaries. We have represented a larger range of 
p to see how these < islands are repeated in the p direction {W can be as small 
as W ~ — 10~^ bordering the limits of our numerical accuracy). 

To see how these structures arise, note that for small K the potential can be 
treated perturbatively [equation (|40|l ] . Then, the wave functions have the form 
^q{x) (X e*«"=(l + Ae+'^ + A*e^'^), with A = i{K/2Tr)'^. The first term is a pure plane 
wave, and the rest first order corrections due to the sinusoidal V{x). Performing the 
Wigner transform of 'i'q{x), we get for the ground state q = (cf. |nSI P- 314]) 

W{x,p)^S{p) + \X\^[S{p-l) + d{p+l)] 

-2|Apcos(2x)<5(p) + 2iAsin(x)[<5(p-i)+(5(p+i)] . (41) 
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The first line gives tlie momentum localisation at p = (unperturbed) and at p = ±1. 
The first term in the second line arises from the interference of e+'^ and e~'^, and 
is responsible for the weak modulation of the height of the p = bell along x. This 
modulation is more intense for larger A (i.e., larger K). The last term accounts for the 
interference between the plane wave centred at p = with those at p = ±1, and can 
produce the negative islands. Therefore, the above simple functional form captures 
most features of the Wigner funtions of figure El 

8. Periodic potentials: quantum dissipative case 

After checking the connection with the classical and Hamiltonian quantum limits, we 
finally include dissipation and temperature in the quantum case. Here examples of 
particles in cosine and ratchet potentials will be discussed. 

8.1. Quantum Brownian motion in a cosine potential 

Transport properties of weakly damped particles in a cosine potential were studied 
by Chen and Lebowitz |^ They started from the system-plus-bath density 

matrix, traced out the bath variables, and performed a perturbative calculation in 
the potential height (followed by a resummation) to get {p). For low forces a free- 
particle like behaviour {p) oc -F/7 was obtained. Increasing F, the wave-vector 
associated to p approaches the first zone boundary (figure |2Jl. There, while Landau- 
Zener tunnelling can bring the particle to the next band, Bragg scattering reduces the 
velocity (cx dE/dq). Eventually, at larger forces, p corresponding to states inside the 
next band become favoured, and the free-particle behaviour is progressively recovered. 

Since high k-^T /h'^ approximations were involved in their calculations, we set 
Dpp — ^Mk^T and D^p = in the master equation Q). Solving it with the continued- 
fraction method we reobtain the effect just described (figure 0] here we work in the 
regime H'y/kBT ~ 0.1). Note that this quantum slowing effect is reduced as 7 increases, 
since the coupling to the bath "broadens" the effective energy levels. This broadening 
makes less relevant the presence of the band gap ("bridges" it) and brings the curves 
progressively closer to the free-damped-particle behaviour (p) oc F/"/. 

With our method we get, in addition, the Wigner function (i.e., the full density 
matrix). The distribution of velocities (|32|l . P(v) = Jda; iy(x,p), when the curves start 
to rise again (i^/7 ^ 4.5), shows two peaks separated by a minimum corresponding 
to the wave vector of the first zone boundary. Loosely speaking the system shows 
coexistence of two quantum "running" solutions. Classical running solutions have 
a wiggling structure along x, as the particle slows down near the potential maxima 
(figure [T]l. The full W{x,p) shows that the quantum running solutions have a nearly 
straight structure (substrate insensitive). We also recognise the familiar crescent 
islands of negative W at large p (there W ~ —0.0003). Around p ^ a. complex 
structure is developed along x ("locked part"), with the maxima deformed and the 
minima becoming slightly negative (W ~ —0.005). These features may refiect some 
interference between the locked and the nearby running solution. However, simply 
adding a plane wave e+'^"^ to the ^'g(a;)|g=o used to analyse the Wigner functions 
of the Hamiltonial problem [equation 141|l ] we have not been able to reproduce such 
structure. Anyway, it is erased when integrating over x to get Pip). 
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Figure 4. Upper panels. Left: (p) vs. F for various dampings at T = 0.5 
and K = 1. Right: marginal distribution of momenta P{p) = f dxW{x,p) for 
7 = 0.01 at -F/7 = 3.5, 4, 4.5, 4.75, and 5. The arrows mark the zone boundaries 
IpI = ■k/K. Bottom panels: two views of W{x,p) for 7 = 0.01 and F/7 = 3.5. 



8.2. Quantum Brownian motion in ratchet potentials 

We now turn our attention to periodic potentials without spatial inversion symmetry. 
This feature combined with out-of-equilibrium conditions can give rise to directional 
motion with net unbiased driving [ratchet effect). This directional motion has been 
invoked to explain the behaviour of molecular motors and, with the emergence of 
nanoscience, opens the way to build mesoscopic devices and engines based on it. The 
theoretical work has concentrated on the classical behaviour and in the high-damping 
regime The few quantum studies also addressed strong system-bath coupling in 
two situations: semiclassical |57| and a extreme quantum case disregarding thermal 
effects [S^. An exception to the large 7 studies is but there the substrate potential 
(and the force) were treated perturbatively. Here, we shall solve the quantum master 
equation Q with the (non-perturbative) continued-fraction method for particles in 
ratchet potentials, taking into account finite damping non- negligible inertia). 

For the ratchet potential (|38|l . the Fourier coefficients of V'{x) are V^i — 
—Vo/2kBT and V^2 = —rVo/2kBT (recall the thermal scaling oiV; section|2}. We use 
r = 0.44, which smooths a small shoulder that this potential exhibits on the "easy" 
side (figure [SJ. As discussed before, the range of index coupling is 2, which gives a 
5-term recurrence. For the stationary response it can be written as [cf. equation 

Qa"C„_2 + + QaC„ + + Qi+C„+2 = (42) 

which can be folded onto a canonical 3-term recurrence by introducing appropriate 
(block) vectors and matrices ( [Appendix A| ). 
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The potential considered is the minimal extension of the cosine potential lacking 
inversion symmetry, so that the average velocities can be different for positive and 
negative forces (vd. figure[Sl). We excite the system with a square- wave force switching 
alternatively between ±i^, and compute the "rectified" current of particles 

l{p)r = 1{V)+F + l(p)-F ■ (43) 

Here (p)±F are the correponding stationary velocities since we consider adiabatic 
conditions {u — + 0). In what follows we briefly investigate the classical limit 
(understudied for finite damping) and then proceed to study quantum corrections. 

8.2.1. Classical case. Figure [S] displays (p)r as a function of T, showing 
the appearance of non-zero rectified velocities (ratchet effect). The rectification is 
optimum at some intermediate temperature (at too low T there is hardly a response, 
and at too high T the potential is irrelevant, so its asymmetry plays no role in either 
case). At 7 = 10 the results coindice with those obtainable from the analytical {p) for 
overdamped classical particles Ch. 11.3] (vd. also [I^IHTIj l 

. 2.r(l-e-^-^/^) 

/o"d:Ee-^(-)/o"dye'^to) - (1 - e-2-^^/^) /J"dxe-'^W/;dye'Afe) 

with 4>{x) = \V{x) — Fx]/T. Inertia, however, broadens the (p)r curves and shifts 
the maxima to higher T (for the lowest 7 = 0.05 the maximum moves to a slightly 
lower T; vd. infra). This broadening is accompanied by a slower decrease of (p)r 
with T. This is consistent with the asymptotic behaviour of {p)r when T — > cx3 
|59|. For overdamped particles it goes as {p)r ~ T~*, while for weak damping it 
decreases only with {p)i- ^ (17/6 < 3). We have checked that our results 

recover these dependences (curves not shown due to the smallness of the asymptotic 
(p), ~ 10-6 - 10-7). 

In our case F < j Fi (the "retrapping" force; section 17?^ . so that at T = the 
attractors are locked solutions. Then, the absolute j{p) vs. T curves start from (p) = 
at T = 0, depart from zero as T is increased and evolve towards "f{p) = F for T 00 
(free damped particle) . In the overdamped case this evolution is nearly a step, whereas 
the slope of the initial raising decreases with decreasing 7 (if F < Fi). To understand 
this, let us assimilate {p) to the escape rate F (modified by some mean free path), so 
that 7(p) ~ 7F. Kramers' theory [SJ shows that 7F increases monotonously with 7. 
Therefore, the "f{p) vs. T curves for finite 7 must lay ordered below the overdamped 
result. On the other hand, the maximum rectification occurs in the temperature region 
where ^{p) transits from to F. This intermediate region is narrower (and the slopes 
steeper) the higher the damping is (figure EJ, while the curves start to raise from zero 
at a lower T. Therefore, increasing the damping the 7(p)r curves will be narrower and 
the maxima shift to lower T, as observed. Eventually, this also allows to explain the 
anomalous curve showing the return (7 = 0.05). Here 7 ~ F, so that at T = the 
system starts close to the zone of bistability. Then the particle can perform incursions 
into the running solution, increasing the initial slope of ^{p)±f but now when the 
damping is lowered, shifting the maximum to a lower T. 



8.2.2. Quantum corrections. After the mandatory exploration of the classical 
limit, we proceed to make the system quantal, by decreasing K. Figure El shows the 
rectified current with K = lb (with 7 bands below the barrier; cf. figure|21) for various 
7, together with the reference classical curves. We have set |F|/7 = 1 is all cases. 
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Figure 5. Top panels: Wigner functions for a ratchet potential in the classical 
limit {K = 1000) with 7 = 0.05 and T = 1. Two forces are used F = -0.1 (left) 
and = 0.1 (right), to show the lack of symmetry in the response (running part 
more developed to the easy side) and the locked regions reflecting the potential 
profile. Bottom left: Rectified classical velocity vs. temperature for various values 
of the damping and force \F\ = 0.05. Right: 7(p) vs. T for the positive force. 



to have roughly the same amount of locked and running components in the solutions. 
At high temperatures {k^T ^ Vq) the classical and quantum curves coincide, as 
expected. However, decreasing T down to fceT ~ Vq, the ratchet effect gets reduced 
or enhanced, with respect to the classical result, depending on the dissipation. Finally, 
at low temperatures {k^T <^ Vq) the rectification is always reduced. 

To understand qualitatively the underlying physics we turn to some semiclassical 
results for various related problems. First, let us consider the reflection/transmission 
coefficient for an asymmetric saw-tooth potential \<o2i § 50, 52]. The results depend 
on the energy of the incident quantum particle. The reflection coefficient for energies 
higher than the barrier, TZ, and the underbarrier transmission coefficient, T, read 

^3-(^2 



n 



(4£:)3 



r 



exp 



dxp{x) 



(45) 



where /i, /2 are the slopes at both sides of the cusp, E the particle energy, and 
p{x) — y^2M[E — V{x)]. Suppose now that we have an asymmetric potential, with 
an easy direction to the right, and we deform it with forces ±i^, computing the 
corresponding reflection coefflcients 7?,+ and TZ- . For energies higher than the barrier, 
the above TZ gives TZ+/TZ- < 1 ( [Appendix F| |. That is, overbarrier wave reflection is 
less intense when the slope of the potential is smaller. Thus overbarrier transmission 
(as thermal hopping) is favoured along the easy direction for energies above the barrier 
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Figure 6. Rectified velocity vs. T with 7 = 0.2, 0.1, 0.05 (using |F|/7 = 1) 
for K = Vo and the reference classical curve (K = 1000). The insets of 7 = 0.2 
and 0.05 show enlargements of the peaks (together with results for K = 20 and 
10). Bottom right panel: Kramers classical rate vs. 7 at various forces (from 
Mel'nikov— Meshkov formula 'BT' ) . The crosses mark the points with F/7 = 1. 



|63|. In contrast, using T for energies below the barrier, we find T^/T- < 1, so 
that > 1. That is, the reflection is more intense in the easy direction or, 

equivalently, tunneUing events are favoured in the hard direction |57j . 

For a distribution of particle energies one would expect a competition between 
these two types of quantum effects. Then, the enhancement of the rectified velocity 
under ±F forcing for the 7 = 0.2 curve follows from the increase of the thermal escape 
with 7 in the intermediate-to-low damping range 7 < 0.5 (figureEJ. At that damping 
more escape events occur and the particles launched over the barrier are subjected to 
the discussed phenomenon of wave reflection, which favours net motion along the easy 
side and hence increases (p)r. Lowering the damping, on the other hand, less escape 
events are produced and wave-reflection becomes progressively dominated by tunnel 
events (favouring the hard direction). At k^T ~ Vq one would expect the tunnel to 
be thermally assisted through the higher bands, which are wider. A simple estimation 
of channel probabilities comes from multiplying a band thermal population (using 
Kramers' rate) by the probability of tunnel (given by the semiclassical transmission 
coefficient T). Here the exponential reduction of T with E compites with the 
exponential increase of T for the upper bands ( [Appendix F| ). At fceT ~ Vq and 
for K = 15 we find that tunnel is indeed favoured through the upper bands. 
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To check that the behaviours found at the maxima are systematic, we increase 
K to K — 20 {'d bands below the barrier) and reduce it to X = 10 (4 bands; here 
the temperatures around the peak are still inside the range of validity of the master 
equation, small h'y /k^T). The results (insets) show clearly the tendency to amplify 
either the enhancement or the reduction when the system becomes more quantal. 

Finally, at low temperatures {kB,T ~ Vb/5) thermal activation becomes ineficient 
for all dampings. Then a reduction of the rectification is expected, due to dominance 
of below-barrier energies and hence of tunnel events (through the lowest bands) . Note 
that the decrease of {p)r relative to the classical results is comparable in the curves 
shown. This suggests that we have similar crossover temperatures Tq in this damping 
range. [Tq, below which quantum effects dominate, is useful as a measure of quantum 
corrections at a given T). This quantity is known exactly in the problem of the escape 
from a "quadratic-plus-cubic" potential V{x) = ^MujIx'^[1 - (2a;/3xo)] Ch. 14] 



To = (fiiJo/27rfcB)(VTT^- a) 



(46) 



In our units Tq (■\/l + 7^/4 — 7/2) //f and in the damping range studied the factor 
in brackets is close to one (^ 0.9, 0.95, 0.975). This gives crossovers To ^ 1/K ^ 0.07 
and hence comparable reductions of the ratchet effect, as observed. This reduction 
of the rectification is the precursor of the current reversals ((p)r < 0) found in these 
systems when tunnelling completely dominates. Although we find small (p)r < in 
some cases, the results cannot be fully trusted. For example, at 7 = 0.2, if = 15 
and T = 0.2 we have h'y/kBT ~ 0.5, which is bordering the limits of validity of the 
quantum master equation employed here. 

9. Periodic potentials: dynamical response 



We conclude with the study of the (non-adiabatic) dynamical response to oscillating 
forces F{t) — AFcos{Lut). Then, the coefficients of the expansion H2()|l of the Wigner 
function are periodic functions of t, so they can be Fourier expanded as follows: 



(47) 



k=l 



To lowest order in the probing field, the static part of the corresponding vectors 
[equation H24|l ] obeys a recurrence relation of the type (vd. [Appendix G| ) 

(48) 



,(0) 



11+ JO) 







The equations for the harmonics have the form (I is the identity matrix) 



iikul + Qa.)ci^^ 



fo 



(49) 



with the "forcing" involving the previous order results. Specifically, 



, with AQq, the part of Qa corresponding to AF, namely 



Vlv+^F 

VT77-AF V2r]+AF 

V2tj^AF VSrj+AF 

V3 77_AF 
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Figure 7. Linear dynamical susceptibility vs. frequency at T = 0.05 and K = 
200 with 7 = 0.01 (~ classical), 0.0003 and 0.0001 (less to more peaked curves). 
Left panel: real part. Right panel: imaginary part (dashed line, 7 = 0.0001 
but halving T = 0.025). Vertical lines: loci of the transition frequencies of the 
associated quadratic-plus-quartic anharmonic oscillator [equation 1511 ]. 



[yd. [Appendix E| for the general Qq/j; the AQq, actuaUy entering in is divided by 
AF, cf. equations l|G.l|l and l|G.4|) ]. Equations I^Hjl and 1^ can be solved sequentially 
with the continued-fraction method. 

FigureElshows the linear susceptibility (fc = 1) computed in this way for a particle 
in a cosine potential as a function of lo. (Here we return to safe ground because 
hj/kBT ~ 10~^-10~^; for related dielectric and Kerr relaxation curves of a quantum 
rotator vd. [Bl] and references therein.) For the largest damping the results correspond 
to the classical limit. There the line-shape x"{^) broadens and extends to frequencies 
lower than the oscillation frequency near the bottom of the potential wells (lu/loq = 1 
in our units). For low damping, this classical spreading of oscillation frequencies has 
an important contribution from the amplitude dependence of the period of oscillation 
in anharmonic potentials 14. 5Q. . This non-dissipative contribution to x"(^) is also 
known in classical spin problems, where the precession frequency depends on Sz due 
to the magnetic anisotropy |65[ 

The quantum regime is approached by decreasing the Kondo parameter a ^ jK 
[equation lO]. Reducing 7 with a fixed K we find that the x(a;) curves develop a 
multipeaked structure. These peaks could be vaguely seen as the mentioned nonlinear 
oscillations becoming quantised. To give content to this statement, we investigate 
perturbatively the effect of the anharmonic terms of a cosine potential (those beyond 
oc a;^) on the energy levels of the harmonic oscillator part. Since the main contribution 
should come from the quartic term, we consider the potential V{x) — ^MloqX^ -\- hx'^ . 
First-order perturbation theory gives the eigenvalues 67, Ch. 6.3] 

^hujo(m+\)+ 36 {h/2Mujof (2™^ -I- 2m + l) . (50) 

Expanding the cosine potential V{x) — — Vb cos(x/a:o), we identify Mujq = Vq/xq and 
b ~ —Vo/{xq4:1) (< 0, soft anharmonicity) . Then we can compute the level spacing 

. \ h (m + 1)1 

A£;„+i,„i = hojQ 1 - — — 2" (51) 

[ SMujqXq 

which has acquired a dependence on the level index m absent in the harmonic case. 
The associated transition frequencies Em+i,m/f>' are plotted in figure [7| [scaled as 
i^m+i,m/'^o = 1 — TT{m + 1)/4:K]. They agree well with the location of the main peaks 
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of the quantum x"('^) curves, supporting their relation with the nonhnearity of the 
potential Furthermore, the contribution from transitions between the higher levels 
should be reduced when they are less populated. Decreasing T, we indeed see that the 
corresponding peaks (located at lower uj) are substantially reduced. In the suggestive 
language of nonlinear oscillations we would say that lowering the temperature the 
large amplitude oscillations (having longer period) become less probable. 

The applicability of the quantum nonlinear oscillator picture is grounded on the 
flattness of the lowest bands at the K considered, which are well approximated by 
constant levels (vd. figure[2Jl. We have computed the semiclassical width of the bands 
in a cosine potential in [Appendix F| An upper bound for the width [equation ljF.9p ] 
reads in dimensionless units Wh = (4/iir) exp [ - ^K{1 - e)] , with £ = E/Vo the band 
mean energy relative to the potential amplitude. For moderate K this bound gives in 
fact very flat bands, specially the deepest ones {E — Vq). This also implies that the 
width of the x"('^) peaks has a negligible contribution from the intrinsic level width, 
being dominated by the damping and thermal broadening. 

10. Discussion 

Traditionally, the continued-fraction method has been employed to solve Fokker- 
Planck equations for few-variable classical systems (translational and rotational 
problems). When compared with numerical simulations of the Langevin equation, 
the method has several shortcomings: (i) it is quite specific of the problem to be 
solved (one needs to recalculate the recurrence coefficients for each potential), (ii) its 
convergence and stability depend on the parameters of the problem (and can be poor 
in some ranges), and (iii) it does not return "trajectories", which in the simulations 
provide helpful insight. Nevertheless, when the method can be used its advantages are 
most valuable: (i) it is free from statistical errors, (ii) is essentially nonperturbative, 
(iii) specially apt to get stationary solutions (static and dynamic), (iv) high efficiency 
(allowing to explore parameter ranges out of the reach of the simulation), and (v) it 
gives the distribution W, which certainly compensates for the lack of trajectories. 

For these reasons, and the lack of quantum Langevin simulations, it was worth to 
develope continued-fraction methods for quantum dissipative systems. This had been 
done for problems of spins in a thermal bath and in quantum nonlinear optics. In this 
article we have discussed in detail the adaptation of the continued-fraction method to 
tackle quantum master equations in phase-space problems, taking advantage of the 
parallels with the classical case provided by the Wigner formalism. We have seen that 
the quantum extension of the continued-fraction method is more problem dependent 
(except for polynomial potentials), due to the necessity of using the a:- recurrence of 
the quantum hierarchies and finding suitable bases. Nevertheless, it inherits from the 
classical methods most of the abovementioned advantages, in particular, the obtaining 
of the Wigner function, from which any observable can be computed. Furthermore, 
the eigenvalue spectrum of the system, although helpful in understanding the physics, 
is not necessarily required. This is advantageous when dealing with non-bounded 
Hamiltonians, continuous spectra, etc. Finally, by changing the appropriate quantum 
parameter, the connection with the classical results is attained in a natural way. 

Note that "continued-fraction" is a quite generic term (like "series expansion"), 
which can be found in many different contexts, and in particular in works dealing with 
quantum dissipative systems (see, e.g., and references therein). However, in those 
cases one usually obtains some continued-fraction expression for a certain quantity 
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(e.g., the linear susceptibility, memory functions, etc.), whereas we have obtained the 
complete solution of the quantum master equation by continued-fraction methods. The 
approach allows to study the interplay of quantum phenomena, nonlinearity, thermal 
fluctuations, and dissipation in several classes of systems. This is important, because 
methods optimised for one of the sides (e.g., nonlinearity or quantal behaviour) tend 
to perform poorly for the others (e.g., fluctuations or dissipation). In addition, the 
visualization of W{x^p) plus the knowledge of the classical phase-space structure 
(orbits, separatrices, attractors, etc.) can provide valuable insight in difficult problems. 
Finally, the main limitations of the approach are those that the starting master 
equations may have. We have considered the celebrated Caldeira-Leggett quantum 
master equation, but the method can be applied to a class of equations of this type, 
and to possible generalisations. 

We have implemented the method in the canonical problem of quantum Brownian 
motion in periodic potentials. We have considered both the cosine and ratchet 
potentials. For the former we have recovered and extended a number of classic results 
and interpreted them on the basis of the Wigner formalism. Results for the dynamics 
under oscillating forcing have also been obtained, illustrating the modification of the 
nonlinear effect of spreading of oscillation frequencies by quantum phenomena. For 
particles in ratchet potentials we have studied, under adiabatic conditions, the effects 
of finite damping on the rectified velocity. Taking into account the competition 
of thermal hopping, overbarrier wave refiection and tunnel events, together with 
semiclassical analytical results, has allowed us to understand in great detail the physics 
of these systems. 
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Appendix A. Solving recurrence relations by continued fractions 

Here we briefly discuss the solution of recurrence relations by continued-fraction 
methods. We first consider 3-term recurrences and then differential recurrences, 
vector-matrix cases, folding of larger-coupling recursions into 3-term ones, and the 
problem of the intial conditions. (For a brief description of the relation of continued 
fractions, series expansions, recursions, and orthogonal polynomials see |69j.) 

Suppose we have the simplest case of a 3-tcrm recurrence relation of the form 

QTC,^l+Q^C^+QtC,+l^-f,, Z = 0,±l,±2,... (A.l) 

with Qi and fi given quantities. To obtain the Ci one inserts in ljA.l(l the ansatze 
[CT ITU] (upper and lower signs for i > and i < 0) 

C^ = S.C^^i + a, (A.2) 

getting the following ladder coefficients Si and shifts 

% , a,.-A±4^. (A.3) 

Qi + Qt Si±i Qi + Si±i 

If the recurrence is finite or the Ci decrease quickly with |i| (e.g., when they are 
coefficients of some judiciosly chosen expansion), we can truncate at some large \i\ = I, 
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by setting S±i = and a±i — 0. Next we iterate downwards in (|A.3|I getting Si and 
Oi, down to i = 0. Then, to get the Ci from ljA.2p . we need Co, which is obtained from 

(Q^^-i + Qo + Q^Si) Co = -{fo + Qoa-1 + Qjai) (A.4) 

combination of equation IjA.ll) at i = and C±i — S±iCo + a±i. Thus, starting 
from Co, we iterate Ci — SiCiz^i + Ui upwards, getting the solution of the recursion. 
For fi = (homogeneous recurrences), we have = and the solution reads 
Cf ^ = SA^i ■ ■ ■ SiC^^\ Note that Si is given in terms of S'i+i in the denominator, 
which can in turn be written as a fraction with 5*^+2 in the denominator, and so on. 
This has the structure of a continued fraction, naming the method 

C= ^ip^. (A.5) 

91 H : 

92 H 

The iterative solution can be cast into explicit form as a series of products of 
continued fractions [7J. The above form, however, has the virtue of being easy to 
implement in a computer. Convergence is checked increasing I ^ 21 ^ ■ ■ ■ and 
repeating the procedure. Solving equation HA.1() in this way requires a computational 
effort of order /, instead of the order P of an arbitrary linear algebra problem. This 
reduction arises from the tridiagonal structure of the matrix associated to (lA.ljl . (For 
analytical inversion of tridiagonal matrices see |72| and references therein.) 

We can solve analogously a differential recurrence of the type 

dcjdt = grc,_i + + q+c,+i + /, . (a.6) 

Laplace transformation converts this equation into [g{s) = J^dte~^*g{t)] 

Q-C,-i + {Q^ - s)C, + Q+C,;+i = + CM] (A.7) 
where g{s) — sg{s) — g{0) has been used. Then, introducing Q'i = Qi — s and 
// — fi + Ci(0), the above equation has the structure of the ordinary recurrence IjA.ip . 

The quantities involved in the above recurrences need not to be scalars, but they 
can be J-vectors (Ci and fi) and the coefficients Qi are J x J-matrices. Then we 
speak of matrix continued fractions. The only change in the above solution is that the 
fraction bars then mean matrix inversion ("from the left", A/B A). 

This is important, because a recurrence relation involving more than three 
coefficients can be "folded" into a 3-term recurrence by introducing vector and matrix 
quantities. Let us show how to do so for the following 5-term recursion 

QrC,-2 + Q^C,-i + Q^C^ + QtC,+i + g++C,+2 = . (A.8) 
Defining the 2-vectors 

c. = ( ) and f , = ( /^^ 

V '^2i+l J \ J2i+1 



and the 2 x 2-matrices 

Q2t Qti \ n- - f Q2i \ (n+ - ( Qtt 



'2i+l Q2i+1 J ' y Q2i+1 J * V '9Jj+l Qti+l 



it can be easily seen that equation IjA.Sp for 2i and 2i + 1 is equivalent to 

Qrc,_i + Q,;c, + Q+c,+i = . (A.9) 

Then, insertion of ~ SiCi^i + ai gives the matrix version of ljA.3|) . Note that the 
Ci and fi can already be vectors and the Qi matrices. Then, the above q, /j, and Qi 
are "block" vectors and matrices. 
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Finally, to iterate upwards — §,;Cj^i + 0^, the "initial condition" Cq is obtained 
from the matrix version of (|A.4p (which we write Ax = b). As in the absence of 
forcing (/^ = 0) the general solution involves an overall multiplicative constant, one 
adds to this system an extra equation to fix it (e.g., a normalisation contidion on the 
coefficients; vd. section (BJ. Alternatively, we can fix arbitrarily one element of Cq and 
rescale the solution at the end. In either case, one has an extra equation which added 
to Ax = b yields the extended (J + 1) x J system 



l.h.s. 



eq. 




/ 



hi 
b-j 



(A.IO) 



\r.h.s. 



This system can be solved by any method appropriate for cases with more equations 
than unknowns (e.g., using QR- or SVD-decomposition [ZSl)) getting Cq. 



Appendix B. Derivation of the transformed evolution operator C 

Here we derive C , the w^^( ■ )w transform of the operator C in the master equation 
with w cx exp[— (77^^/2 + $)] [equation ifTUIl ]. We split £ as £ = £ir + £kin + £v, 
with £ii-, £kin, and £v given by (jSJ. For the "product" of operators we use that the 
transformation of the product is equal to the product of the transformations: 

TB^AB, A^^CA)"" (B.l) 

which follows by sandwiching "identities" ww~^ — 1 between the factors. Besides, 
the multiplication by functions of {x,p) commutes with w~^{-)w, so that only the 
differential terms need to be transformed. From these two properties we see that only 
dx and dp need to be calculated (and then raised to some power or multiplied). 

Appendix B.l. Calculation of dx and dp 

Calculation of dx. In order to transform dx, we apply w^^dxW, to an arbitrary function 
f{x,p). The p dependent part of w cancels out and we obtain 

d^f^e^i-dx'^e-^f + c-^dxf) , ^ d^ ^ dx ~ ^' (B.2) 

which has the form of a "displaced" differential operator. 

Calculation of dp. Applying w^^dpW to a function /, we similarly get dpi 

^/ = e''P'/2(-r,pe-"P'/2/ + e-''P'/'ap/) , ^ Wp^dp-T^p. (B.3) 

This dp can be expressed in terms of the following creation and annihilation operators 
h=dp + \p, h+ = -dp + \p. (B.4) 

Thus, using the shifted rj parameters ri± — r/ ^ ^ [equation H14|l ] we arrive at 

d;=-{rj-b+ +7]+b) . (B.5) 

For rj = 1/2 (the choice in the classical case), r/^ — 1 and = 0, so that dp = — 5+. 
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Appendix B.2. Calculation of £[,■ ~ w^^Cn-w 

Using the product property AB ~ AB, we have Ci^ = ^rdpiv + dp) ■ Then, replacing 
dp by its "bar" form (|B.5ll and writing p = b +6+ [from (|B.4|) ]. we find 

-7y^£ii. = ry+(f - + ?7_(1 - r]^)b+b+ + 2{rj - r]^ri+)b+b + ?/+(f - 7]+)bb . 

Here we have also used + ry+ = 2ry and the commutation rule [b , 6+] = 1 to get a 
normally ordered form. If we introduce some compact notation for the coefficients 

7d = 2-fT {ri - i]-r]+) 7_ jtV-C^-V-) 
7+- = 7T ??+(!- J?-) 7+ = 7T?7+(l-77+) 
we can finally write 

Ar = -{-/+-+ 7-b+b+ +-f,ib+b +j+bb) . (B.6) 
(The choice rj = 1/2 gives the self-adjoint form = — 7t |14j.) 

Appendix B.3. Calculation o/ £kin = w^^'CkinW 



Using equations ljB.2|l and ljB.5|l for dx and dp, as well as p = 6 + 5+ and the definition 
H15|l of d± we get for the kinetic part = — (p — D^p'O^ dx- 

£ki„ - - $') - - $') . 

The coefficients of b and 6+ are Z)^ and also defined in p5|) . and we have 

£ki„ = -(6i^+ + 6+^") . (B.7) 

Appendix B.4- Calculation of Cy = w^^CyW 



To obtain Cy = X;s>o '^^'^''^^^'''^^^^p^'^^^ equation ©], we use the product 
property IjB.ip to raise 9p to 2s + 1, getting 

= _ ^ k(^) 1/(2^+1) + 77+5)'^+' . (B.8) 



s>0 



To get later the matrix elements of Cy between Hermite functions it is convenient to 
cast (|B.8|I in a normally ordered form. To this end we use Pathak's theorems [7^ in 
the following generalised form 

/ I \2s+l s ./J- n2(s~(7) + 1. , , , 

(«++«) ^ V •(" +") • (J. n^ 

(2s + 1)! ^ [2(s-g) + l]! q\ ■ ^ ■ ' 

Here :(a+ +a)™: EfcU CfcH^"^ with 4" = m!/[fc!(m- fc)!] ( :/: is the 
result of moving the a"'' to the left as if they were scalars) . The generalisation resides 
in letting a and a"*" have a constant but different from one commutator. In our case 

a=r]+b, a+=f]^b+, [a,a+]^cr (B.IO) 

follows from the commutator [b , 6+] = 1 with a — ri^rj^. Note that a and are not 
the adjoint of each other, but the notation is symbolic. 

With help from equation ||R9)| and k'^"^ = (-l)'*A^^/(2s + 1)! [equation ©] we 
can write (|B.8|I as follows 

oo s 
s=0 9=0 
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where we have introduced the notations 

[2(, -,) + !]! q\ ^ ' X-.[a+a).. 

For a given s different powers of X appear. We take a common factor and add 
coefficients, e.g., Aqq + An + A22 + ■ ■ ■ with X, Aio + A21 + A32 + ■ • •, with X^, etc. 



00 00 



s=0 9=0 £=0 s=i i=0 q=0 

For s = (7 + in As^q we write _ ^2qy(2^+i) ^^^^^ introduce the operator A 

(-A/2)'^ (-l)^^dB_T/(2£+l) A — c)2 m 19^ 

A+f.? ^1 (2£+ 1)! ' ^-'^-^dBC'x- (^-12) 

Now the q and ^ dependences factorise, so we can sum the series of q in Ei = 
^^Q^q+£.q, getting an exponential. Then, recaffing the definition of k''*^ we finally 
obtain the normally ordered form for Cv we were looking for 

oc 

£, = -Y,Er.{a+ + ay'^': , = «;(^) exp(-A/2) W^^+D . (B.13) 

i=0 

The operator exp("A/2) does not act on what could appear after V{x), since A was 
introduced as some abreviated notation for the high-order derivatives of the potential. 

Appendix C. The different contributions to dtCn = J2m QnmCm 

Here we calculate the matrix elements Qnm = Jdpipn^ijjm, with the transformed 
operator £ = £ir + £kin + ^-v obtained in [Appendix B| We split in parts the 
calculation by introducing the corresponding decomposition: Qnm — ("■I'^irl'Ti), 
Qnm = ('T-l'CkinIm-), and Q^„„ = (ri|£v|"i), where we have used the custom "bra"-"ket" 
notation. The calculation is simplified because we have all operators L expressed as 
normally ordered forms of h and 6+ . Then we can take advantage of the ladder actions 
h^\ri) = VrT+l |n + 1), b |n) = ^/n \n -— 1), and the number property b^b \n) = n\n). 

Appendix C.l. Calculation of J2,nQnmC^m 

From ljB.6|l for Ci^ and the orthonormality of the |m), we find for QJ^m = {n\Cir\m) 



-Qnm = 7- ~ '^)nSn-2,rn + (njd + 7+-) ^nni + 7+ V (^^ +!)("- + 2) 5n+2,m ■ 

Then, we introduce the notation 7„ = n 7d + 7h for the coefficient of the diagonal 

term [cf. equation multiply by Cm and sum over m, getting 



fnmCm = 7- V - l)n C„_2 + 7n C„ + 7+ \/in+^){n+2) Cn+2 ■ (C.l) 



Appendix C.2. Calculation of ^^_^_^Q^^C, 



Using equation (|B.7|I for £kin and the "ladder" action of b and 6+, we have 

Qnm = ~{n\bD+ + b+D-\m) = - (xATTT D+6n+l.m + ^ D-6n^l,m) 

where we have used the Kronecker delta to interchange n and m=p 1. Then, multiplying 
by Cm and summing over m, one obtains 



^kin 

Inm^rn 



Cm=-{V^ D-Cn-1 + D+Cn+l) ■ (C.2) 
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Appendix C.3. Calculation ofJ^mQmn^ri 



This will take longer. First, the operator !(a++a)^^+^! appearing in equation (jB.13|l 
for Cv can be put in explicit form by using the binomial formula 

21+1 

■.{a+ + af^': = cf+\a+f a' . (C.3) 

The a and a"'" are proportional to b and 6"*" [equation (|B.10|l ]. so the contribution of 
each summand to {n\Cv\m) involves {n\{b'^)'^^'^^~''b'^\m) cx (n — {2£ + 1) + k\m — k). 
From orthonormality only one term of the sum ljC.3|) will contribute, namely 2k = 
{m — n) + 2£ + 1, while m — n will be odd, m — (2s + 1) with s > (reflecting the 
odd powers of dp in C^). Then, in terms of s, the index k is restricted to 

r £ - s , for m = n-{2s + 1) 

\ ^ + s + 1 , for m = n + (2s + 1) ' ^ ' 

Relating the results for m = n ^ (2s + 1). The two cases can be easily related. 
First, the combinatorial coefficients are equal c^+^^i — '^i-gi because c™ — c™_^. 
Second, the r\ factors, coming from the proportionality of the a, a"'', to b,b'^, are 
simply obtained from each other by exchanging pluses and minuses (cr = rj^rj^): 

fc = ^ + s+l - r^i-'r^i+-'+' = a'-'7jl'+' ^^■^> 

Third and last, the matrix elements are also relatable. For m = n — (2s + 1), we have 
{n\{b+Y+'+'^b'^-'\n - (2s + 1)), while for m = n + (2s + 1) 

{n\{b+Y~''h'^+'+^\n + (2s + 1)> = {m\{b+Y+'+^b'^-'\m - (2s + 1)) 

which is formally identical to the former by replacing n n' = m and m m' = n. 

Thus, from now on we consider only the case m = n — (2s + 1) (factor cr^~'^r]'^^^). 
The result for m = n + (2s + 1), will be readily obtained from this by replacing 

n' = m , m' = n' - (2s + 1) and r^^s+i ^ ^2s+i ^ q-^ 
From all these considerations, we can write the matrix elements of !(a^ + a Y^'^^'. as 
{n\ :(a++ a \m) = 7^^J+^cf_+' a'-' {n\ib+Y+-'+'b'-'\m) 

where we omit a Kroneckcr delta ensuring m = n — (2s + 1). The rest of the calculation 
consists of multiplying this by Eg and sum over all £ to get (n|£v|"^) [equation (jB.lSp ]. 



Restrictions on the sum on £. The contribution of the infinite sum (|B.13p to 
(n|£v|"^) is fortunately cut-off. First, the lower index in c^J^Y should be positive and, 
introducing the shifted index £' = £ — s, we obtain (we rename ^' ^ ^ at the end) 



(n|£v|m) - -rj^J+'J2E^+^ cf '+'^+' a' {n\{b+Y^^+'^^V\m) . (C.7) 

Here we have already used that the sum is also cut-off from above at ^,nax — 'ti, since 
b^\m} = 0, for £ > m. On the other hand, for Eg^s [equation (|B.13|) ] we have 

P -A/2 (2s +1)! 2 si2\<' (~l)''^dB 1^(23+1) /p ON 

^'+^-' [2(^ + s) + l]!("^'^i^^-) (2s -Kl)!^ ■ ^^-^^ 



Continued- fraction methods for the Wigner Caldeira-Leggett master equation 29 

Combining this with c^*^^^*-*^^, gathering (A^gf?^)^ with to form the operator A 
[equation ljB.12|) ]. and recaUing the definition © of k!''^\ we obtain 

^-e+V(-)e-A/^f pl+y+^l! ("|(^")^""^^"^^» ^ ^'""^^ 

Note that the quotient of factorials equals l/(2s + 2)^, with the Pochhammer symbol 
defined by ,3L Ch. 13.5] 

{a)i = a(a+ 1) •• - (0 + ^- 1) = (a + £ - l)!/(a - 1)! . (C.9) 

Computation of the matrix element. To conclude, we only need the matrix element 
(n|(6+)(2''+i)+^5^|m) in We first pass 2s + 1 times the action of b+ to |n) by 

taking its adjoint h and then we use repeteadly its downward ladder action, getting 

h'^'+^\n) = ^/n{n-l)---{n-2s) |n-(2s + l)) = ^/n\/m\ \m) (C.IO) 
with m = n — (2s + 1). Then, in terms of the Pochhammer symbol (|C.9|) . we can write 

= y^ny^ {m\{b+Yb'^\m) = V'^ny^ {-l)\-m)i . 
Now, taking into account that the factor (—1)^ cancels that of (— A)^ and using the 
modified coefficient Kn'' — n^"'^ n\ / m\ [equation we arrive at 

(-m), A^ 



r m 

l^ (2s + 2), £! 



y(2s+i) 



Comparing with the series definition of the confluent hypergeometric (Kummer) 
function iFi(a, c;z) |^ Ch. 13.6] we get 

(n|£v|m)=-e+i4^)[e-^/S^i(-m,2s + 2;A)]l^(2^+i) ,F^{a,c,z) ^Y.Ty\ ■ 

The operator in front of V(x) is simply P^'" [equation (|16|l ]. so that we can finally 
write 

Q«m = (^Pv|m) = -P,^r^^''+'Ha;) , m = n - (2s + 1) . (C.ll) 

Note that \F\{—m, 2s+2 ; z) is a polynomial because the first negative index cuts the 
upper Pochhammer symbol. Indeed the "integer-plus-one" second argument (2s+l) + l 
tells us that it is an associated Laguerre polynomial Ch. 13] 



mz)='-^^^F^{-l,k+\-z). (C.12) 
In our case I = m and /j = ra — m = 2s + l,so that \F\{—m, 2s -I- 2 ; z) = i^*+-'^(z)/c™. 



Final expressions for J^mQnm^rn- Once we have computed Qnm^ ^^^^ matrix 
elements of £v, we multiply by Cm and sum over m (i.e., over s), to get 

7n<n s>0 

Note that this s-sum is restricted because n— (2s -I- 1) > 0. This yields the upper limit 
Smax — [{n — l)/2], with [a] the integer part of a. The contribution from m > n is 
easily obtained recalling the transformation rules (|C.6|I (here the sum is not restricted) 

m>n s>0 
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with r*'+ given by Hlt)|) . We have enclosed into square brackets the action of the 
operators F^'* to recall that they only act on the x-dependence of the potential. 

For 77 = 1/2, the choice in the classical case, rj- = 1, 77+ = 0, and cr = 0. 
Then, r^'^(o- = 0) = and r^'+(tT = 0) = 0, so that only the C„_(2s+i) terms 
survive. The chain of equations then adquires an unbalanced structure (like in the 
old hierarchies |38l I40L HTj l. which may result in a poor stability when the quantum 
terms are important. In constrast, 77 ~ gives \ri±\ ~ 1/2, and hence the weight of 
the terms at both sides of n above and below the diagonals in the matrices Qq^) 
is similar, improving the stability when going into the deep quantum regime. 



Appendix D. The auxiliary integrals 1^ and Kn^ 

For special r/ and the integrals and Kn'' [equation appearing in the 

expressions (|30|) for the observables, are easily done. For ry = 1/2, we simply have 
k!^^ = 5nfl, kI^^ = and K^n^ = 5„^o + ^/25n,2■ If $ = (i.e., e = 0), we 

get la = la^afi- Then the first moments reduce to 1 = IqCq (normalisation) and to 
{p) = IqCi and (p^) = Io{Cq+\/2C2)- In the general case, one can derive recurrences 
relating the /C^^'' with lower order (in £) ones. Thus, only Kn^ is needed, which can 
be found analytically. Concerning the la, we mostly use £ = 0, and hence la — loSa.o- 
Anyway, for periodic potentials one can derive recurrence relations for them which 
can also be solved by continued fractions. 

Appendix D.l. The integrals Kn^ (arbitrary potentials) 

To get the recurrence for the kIP we start from their definition (|29ll with index 
£ + 1, express the last p as p — 6+6+ [equation (|B.4|) ]. and then use pipn = 
y/n^Jn-i + y/n + 1 V'n+ij getting 

= V^Ki'l^ + V^Ki%, . (D.l) 

This recurrence can be used to get the Kn from the Kn ■ 

We can get analytical expressions for the first few Kn with help from the 
tabulated integral 75, equation (7.374-4)] 

'dxe-^' H^{x)H2n+m{ax) + (^a)"^ (a^ _ 1)" . (D.2) 

For Kn^ — "^Q j e~^P 4'n, we use this result with m — 0, obtaining 



Kiti = , Ki°J = V2 A« , A = -,+ /2,_ . (D.3) 

Although Ki^^ and K^^^ can now be obtained from ljD.l|) , they can also be done with 
help from equation (|D.2|I . getting: i^2^'' — 0, i^2n+i — 0, and 

,(1) /TTT— TTT , /77 



n'. ^-^'^ "' n\ j^-^i^ \ri- 
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Appendix D.2. Recursions for the integrals la (periodic potentials) 

For a general periodic potential we can derive a recurrence relation for the integrals 
la — /o^'^dx e-xjp[~eV{x)]ua{x) [equation it^ ]. We integrate by parts the expression 
for ia la and use V = e'"^ with Ua — e'"^/\/27r, getting 

ia la ^ e^VpIa+p ■ (D.5) 

For a finite number of harmonics (V^ = 0, if > B), the recurrence (|D.5p has a finite 
coupling range. For instance, for the ratchet potential H38II . we have V^i = —Vo/2kBT 
and Vj^2 ^ —f'Vo/2kBT and the above recurrence reduces to 

iaIa = -eT[(Ia-l+Ia+l)+r {Ia-2+Ia+2)] , £t = £ Vb/2fcBT . (D.6) 

With Iq as input (obtained numerically, e.g., by Simpson rule 37, App. 2]), this 
recursion can be solved by continued fractions | |Appendix A| ). 

Appendix E. The matrices Qq/j for general periodic potentials 

Here we write expliticly the matrices Q^/j for arbitrary periodic potentials. We shall 
give the results for a slightly generalised basis, obtained by shifting the origin of 
momenta by an amount p [cf. equation (|26|l ] 

W{x,p) = VC^ Ua{x) 4>n{p - p) ■ (E.l) 



This momentum shift is convenient to "catch" solutions centred far from zero p, where 
the ordinary basis would need many ipn{p) to reconstruct the distribution. (A useful 
choice is, when generating a curve, to set p = (p) of the precedent point.) The 
momentum shift is handled as a change of variable p pp — p — p. Then, dp — dp^ , 
but the term p dx makes the force to enter in the combination (scaled units) 

Fp = F--iTp. (E.2) 

As mentioned in section Q we extract the force F from V{x) and set the auxiliary 
potential $ proportional to the periodic part $ = eV . 

The Qa/3 are obtained by inserting the formulae for D^^, and \r'^^V^'^'^^^'^]ai3 
[equation (|35|l ]. into the general matrices of section |31 Instead of giving a single 
expression, we write separately the matrices for [3 = a and /3 = a ± g. The central 
matrix Qqq is determined by the terms involving 5a^3^ and reads 

^ 7o+pick \/T(i(i+a — V 1 -2 7+ 

^/\{\d-a — ri^Fp) ^i+p\a V2{i d+a-i]+Fp) V2~3-f+ 
\/l • 2 7- V2{id-a-r]^Fp) 72+piQ! V3{id+a-r]+Fp) V3 • 4 7+ 
-\/2 • 3 7_ \/3{id-a — ri^Fp) 73+piQ! -\/4(i d+a — ry+i^p) 

V3 • 4 7_ V4(id-a-?/--Fp) "fi+pia 

V 

with 7-1- and 7„ given by H13|) . r/± = and d± ^ 1 + ri±Dxp. Here the periodic 

potential does not appear; only F enters in the matrices Qaa- Besides, they do depend 
on the index a (cf. below). 
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The matrices Qa.aiq are determined by the terms involving V^_p in equation ()35|l 
and read 



,a±q 



= V- 



T9 



-pe r^+Gl 







-^/le 
pe 71+Gl 



-\/2£ 







„5 r<2 



7y_G^-V2e 






-~pe 
T]-G'i~V3e 




V 















G2 

3G1 



fl+Gl-V^e 



n 



-pe 7]+Gg-> 
_G^-V5e -pe 



be 



V 



where G* is evaluated at z = — crq A^g, with a — rj^rjj^. These matrices enjoy 
properties opposite to those of Q^q • The substrate potential does appear (multiplying 
all elements in front of the matrices) while the force is absent. Besides, the matrices 
Qa,Q±(j are independent of the index a, which permits to write them simply as Q^'?. 
The only technical complication with respect to the common choice r] ~ 1/2 is that 
we need to compute the Kummer functions iFi[a, c ; z) included in G^, which in our 
case are simple polynomials [equation ljC.12p ] . Note that to compute the observables 
we must eventually undo the momentum shift Pp ^ p = Pp + p. 



Appendix F. Semiclassical analytical results 

In this appendix we derive some semiclassical formulae discussed in the main text. 



Appendix F.l. Reflection and transmission for a saw-tooth harrier 

In section [8.21 we invoked the current through a saw-tooth barrier deformed by ±Fx 
(figure 1^3). We introduced the reflexion coefficients 7?.+ and TZ- for ±F. One has 
to analyse separately the cases of incident particles with energies above or below the 
barrier. In the first case, the reflexion coefficient by a potential with a cusp (with 
slopes /i and /2) is given by the first equation H45|) '62' § 52]. Owing to the right 
slope in figure lFll is infinity, we make it finite and finally take the vertical limit, getting 

n+ _Vo-FL 



Vo + FL 



< 1 



(F.l) 



Next we consider the case of tunnelling through the barrier. The semiclassical 
transmission coefficient is also written in (|45ll . Defining A = \ jl^dxp{x)\, we calculate 
and A_ for positive and negative forces. After simple integrals, we find 



A^ 



Vo + FL fVo-FL-E 



3/2 



(F.2) 



Vo-FL \ Vo-E 

with E the particle energy (cf. reference @H1)- For moderate forces, F < 0.618 Vq/L, 
we get for the transmission coefficients T = exp[— (2/fi,)A] 

A+/A_ > 1 r+/r_ < 1 . (F.3) 

Thus, for E below the barrier TZ+ /TZ- > 1, while 71+ /TZ- < 1 for E above the barrier. 
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Figure Fl. Left panel: profile of the saw-tooth potential barrier for non- 
deformed (a) and deformed by —Fx (b) and by -\-Fx (c). Right panel: width 
of the energy bands W\, in a cosine potential (exact and various approximations). 



Appendix F. 2. Transmision coefficient and energy hands for a cosine potential 

In order to get the semiclassical transmision coefficient through a potential V{x) = 
— Vb cos(x/a;o), we need to do the integral 



dx y^2M[E + Vo cos(x/a;o)] 



(F.4) 



with turning points obeying cos(a/a;o) — cos(6/xo) — —E/Vq. Due to E + 
Vb cos(a;/a;o) < below the barrier, we change the sign inside the square root 
extracting an i which is absorbed by the modulus. Simple tranformations bring the 
integral into a tabulated form |7S1 equation (2.576)], and we obtain 



A = ASo [2E(s)- (l + e)K(s)] 



E/Vo , 



(F.5) 



where 5*0 = \/ MVqXq while K(s) and E(s) the first and second complete elliptic 
integrals [33 Ch. 5.8]. Insertion in T = exp[~{2/h)A] gives the transmision coefficient 
through a cosine potential. This T also gives the semiclassical energy bands '76' § 55] 

Eg^ Eo± {huJo/TT)VT cos(27ra;og) . (F.6) 

with loq the oscillation frequency in the wells. The width of the bands is therefore 

W^ = 2{hu;o/7T)exp{^4{So/h)[2E{s)-{l+e)K{s)]} , e^E/Vo. (F.7) 

The parameter s goes from 1 when e = — 1 (potential bottom) to at the barrier top 
e = 1. For s ^ 1 we can expand the elliptic integrals, getting 

^ 2{hwo/7r) exp [ - (So/h) (s - (1 + £){l + In [^(1 + e)] } 

+ (l + e)2{l + iln[3L(l + e)]})] . (F.8) 

At the potential bottom e = — 1 this gives Wb — <3xp[~8 (So/h)]. The band width 
increases exponentially with the energy. Near the top, e ~ 1, one has s ~ 0, and we 
can use the series expansion of the elliptic integrals, obtaining 

~ 2{hiuo/n) exp [ - niSo/H) (1 - e)] . (F.9) 
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Figure [FTI shows that equations Hi^'.8() and (|i^'.9|l approximate well the exact results in 
their respective ranges (the later also serves as un upper bound). If we disregard the 
terms with (1 +e)^ in ljF.8|) . the result happens to fit remarkably Wb in all the range. 



Appendix G. Perturbative solution for periodic forcing 

In this last appendix we shall solve the following generic differential equation 
dr 

r- + Qc^f{t)Vc, /(t) = /+e+-* + /_e— ' (G.l) 

perturbatively in the forcing /. (For cosine forcing f+ = f-= //2.) Comparing with 
section |31 we see that c corresponds to the vectors c„, Q to the static part of the 
matrices Qa/j, and the perturbation V to AQo,q, the part of involving /S.F. 
A periodic solution can be Fourier expanded as 

oo 

c = + {fl C^''')c+''^''t + f^ c^k)^-iu^t^ ^ uj^^kuj (G.2) 
fc=i 

with C'-*'-' the unperturbed response, while C*^*^-' ^ c'^^^* for complex c. Separating 
terms oscillating with e+'"''* and e"'"*"*, the left-hand side of (jG.l|) reads 

1 oo 



k=l 



+ ^ [ - i(LjfcT)C('=) + gc(fc)]e-i-^-* . 

k=l 

To obtain the right-hand side, we use ± w = kuj ± w = w^ii, redefine the indices 
(keeping the same names) to get all the oscillating factors at iwfef, and introduce the 
definition C^") = C^o), arriving at 

oo 

/(t)Vc = /+/_(VC(i' -t-VC(i)) ^ ^/^[vcC^-i) 4- (/+/_) VC(''^+i)]e+'"^* 

fc=i 

oo 

fe=l 

Equating terms with the same oscillating factors at both hand sides {uniqueness of 
the Fourier expansion), we obtain 

QC(«) = V[0 + (/+/-) (CW+C(-i))] 

i(a>fcT)CW + QCW = V[C('=-i) + (/+/-) CC^+i)] 

The CC^) are included as C(-'=) = C**^) with = -Wfc. 

The contribution VC'^'''^^^ comes from products of "rotating" terms e='='"* x e='="^'=*. 
For instance, the oscillating factor c+"^*-i* in C''^*^"^^ when multiplied by the part e+"^* 
of the field raises the harmonic from fc— 1 to fc. The term VC^''+^-' comes from products 
of "counter-rotating" terms e^'"* x 6^"^*=*. For example, the product of the oscillating 
factor c+^'^k+it g^jj^j ^j^g g-it^t pj^j-t; of the field lowers the order from fc + 1 to fc. This 
contribution (the "contamination" from higher harmonics) has an order higher in / 
than the contribution of C*-'^"^-', which is therefore the leading one. Indeed, expanding 
the C'^''^ in powers of / we get the leading terms at each harmonic 

Q C(°) = , i{LOkT)C'^''^ + Q CC^) = V CC^-i) . (G.4) 
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The fc = equation gives the static response [corresponding to ()48|l ] , k — 1 the hnear 
dynamical response, etc. [corresponding to equation (|49|l ]. 
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